#setwd("filepath")
path<-getwd()
library(foreign)
library(rjags)
set.seed(666)

### ECRT data
ecrt.dat<-read.table(file="ecrt.csv", sep=",", header=T)
ecrt.dat<-ecrt.dat[ecrt.dat$year>1945,]

### Democracy data
dem.dat<-read.table(file="uds_summary.csv", sep=",", header=T)
dem.dat$cowcode[dem.dat$cowcode==315]<-316
dem.dat<-dem.dat[dem.dat$year>1945,]

new.dat<-merge(ecrt.dat,dem.dat,by.x=c("ccode","year"),by.y=c("cowcode","year"))

new.dat$new.panel<-NA
for (i in 6:dim(new.dat)[1]){
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-1], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-2], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-3], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-4], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-5], 1, 0)
}
new.dat$new.panel[1:5]<-0

post.draws<-matrix(nrow=dim(new.dat)[1], ncol=10)
for (i in 1:dim(new.dat)[1]){
	for (j in 1:10){
		post.draws[i,j]<-rnorm(1, mean=new.dat$mean[i], sd=new.dat$sd[i])
	}
}

difs<-matrix(nrow=dim(new.dat)[1],ncol=50)
for (i in 6:dim(post.draws)[1]){
	for (j in 1:10){
		difs[i,j]<-post.draws[i,j]-post.draws[i-1,j]
		difs[i,j+10]<-post.draws[i,j]-post.draws[i-2,j]
		difs[i,j+20]<-post.draws[i,j]-post.draws[i-3,j]
		difs[i,j+30]<-post.draws[i,j]-post.draws[i-4,j]
		difs[i,j+40]<-post.draws[i,j]-post.draws[i-5,j]
	}
}

new.dat<-cbind(new.dat,post.draws,difs)
names(new.dat)[14:73]<-c("draw.1","draw.2","draw.3","draw.4","draw.5",
"draw.6","draw.7","draw.8","draw.9","draw.10",
"draw.1.dif","draw.2.dif","draw.3.dif","draw.4.dif","draw.5.dif",
"draw.6.dif","draw.7.dif","draw.8.dif","draw.9.dif","draw.10.dif",
"draw.1.dif.2","draw.2.dif.2","draw.3.dif.2","draw.4.dif.2","draw.5.dif.2",
"draw.6.dif.2","draw.7.dif.2","draw.8.dif.2","draw.9.dif.2","draw.10.dif.2",
"draw.1.dif.3","draw.2.dif.3","draw.3.dif.3","draw.4.dif.3","draw.5.dif.3",
"draw.6.dif.3","draw.7.dif.3","draw.8.dif.3","draw.9.dif.3","draw.10.dif.3",
"draw.1.dif.4","draw.2.dif.4","draw.3.dif.4","draw.4.dif.4","draw.5.dif.4",
"draw.6.dif.4","draw.7.dif.4","draw.8.dif.4","draw.9.dif.4","draw.10.dif.4",
"draw.1.dif.5","draw.2.dif.5","draw.3.dif.5","draw.4.dif.5","draw.5.dif.5",
"draw.6.dif.5","draw.7.dif.5","draw.8.dif.5","draw.9.dif.5","draw.10.dif.5")

for (i in 24:73){
	new.dat[,i][new.dat$new.panel==1]<-NA
	}

new.dat<-new.dat[new.dat$year>1949,]

civwar.data<-read.csv(file="ucdp_onset2012.csv", sep=",", header=T)
cw.dat<-data.frame(cbind(civwar.data$year, civwar.data$gwno, civwar.data$incidencev412))
colnames(cw.dat)<-c("year", "gwno", "conflict")
cw.dat<-cw.dat[order(cw.dat$gwno, cw.dat$year),]
cw.dat$gwno[cw.dat$gwno==315]<-316

newer.dat<-merge(new.dat,cw.dat,by.x=c("ccode","year"),by.y=c("gwno","year"),all.x=T)

legal.dat<-read.dta(file="mitchell_jpr.dta")
legal.dat<-data.frame(cbind(legal.dat$cowcode, legal.dat$year, legal.dat$common))
colnames(legal.dat)<-c("ccode","year","common")

newest.dat<-merge(newer.dat,legal.dat,by=c("ccode","year"),all.x=T)

## fill in missing legal system data
newest.dat$common<-ifelse(is.na(newest.dat$common), 0, newest.dat$common)
common<-unique(newest.dat$ccode[newest.dat$common==1])
for (i in 1:length(common)){
	newest.dat$common<-ifelse(newest.dat$ccode==common[i], 1, newest.dat$common)
}

save(newest.dat,file="ecrt_data_NAs.rda")

dat<-newest.dat[,-6]
dat<-na.omit(dat)
dim(dat)

### dropped due to lack of data: Liechtenstein, Monaco, Andorra, San Marino, Slovakia, Montenegro, Kosovo
### Belarus included though not CE member

## Create id for indexing random effects in JAGS
dat$end<-ifelse(dat$year==1998,1,0)
dat$end<-ifelse(dat$ecrt==1,1,dat$end)
dat$end[dat$ccode==265&dat$year==1989]<-1
dat$id<-rev(cumsum(rev(dat$end)))
dat<-dat[order(dat$id, dat$year),]

## Create time counter
helper<-rep(1,dim(dat)[1])
dat$t<-as.vector(unlist(tapply(helper,dat$id,cumsum)))
#dat$t2<-dat$t^2
#dat$t3<-dat$t^3

common.l2<-as.numeric(tapply(dat$common, dat$id, mean))

the.hi.model.list<-list()
for (i in 1:10){

	N<-nrow(dat)
	id<-dat$id
	ecrt<-dat$ecrt
	dem<-dat[,i+12]
	conflict<-dat$conflict
	common<-common.l2
	t<-dat$t
	t2<-dat$t2
#	t3<-dat$t3

	the.data<-list("N"=N, "id"=id, "ecrt"=ecrt, "dem"=dem, "conflict"=conflict, "common"=common, "t"=t)

	the.inits.1<-list(b=c(0, 0, 0, 0), a=rep(0,41), tau=1, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits.2<-list(b=c(1, -1, -1, 0.01),a=rep(-1,41), tau=5, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits<-list(the.inits.1,the.inits.2)

	the.hi.model<-jags.model("ecrt_hi_logit.bug", data=the.data, n.chains=2, n.adapt=0, inits=the.inits)
	update(the.hi.model,90000)
	the.hi.model.list[[i]]<-coda.samples(the.hi.model, c("a", "b", "g"), 10000, thin=1)
}

lapply(the.hi.model.list,summary)
lapply(the.hi.model.list,gelman.diag)

save(the.hi.model.list,file="ecrt_demsq_models.rda")

the.hi.model.dif.list<-list()
for (i in 1:50){

	N<-nrow(dat)
	id<-dat$id
	ecrt<-dat$ecrt
	dif<-dat[,i+22]
	conflict<-dat$conflict
	common<-common.l2
	t<-dat$t
	t2<-dat$t2
#	t3<-dat$t3

	the.data<-list("N"=N, "id"=id, "ecrt"=ecrt, "dif"=dif, "conflict"=conflict, "common"=common, "t"=t)

	the.inits.1<-list(b=c(0, 0, 0), a=rep(0,41), tau=1, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits.2<-list(b=c(1, -1, 0.01),a=rep(-1,41), tau=5, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits<-list(the.inits.1,the.inits.2)

	the.hi.model<-jags.model("ecrt_hi_logit_dif.bug", data=the.data, n.chains=2, n.adapt=0, inits=the.inits)
	update(the.hi.model,90000)
	the.hi.model.dif.list[[i]]<-coda.samples(the.hi.model, c("a", "b", "g"), 10000, thin=1)
}
 
save(the.hi.model.dif.list,file="ecrt_dif_models.rda")

lapply(the.hi.model.dif.list,summary)
lapply(the.hi.model.dif.list,gelman.diag)



